!***************************************************************************
! This module contains the function prob_wage used in E_likegaussj
!***************************************************************************
! This version verifies that at convergence the conditional probability of assignment to A's job levels ~ 1   


MODULE mod_wage


USE mod_types
USE mod_parameters

INCLUDE 'link_fnl_static.h'
USE ANORDF_INT

IMPLICIT NONE


CONTAINS


!********************************************************************************
!Probability of the observed wage with kernel exp(- 0.5 *((X - MEAN)/SD)**2)
!********************************************************************************
FUNCTION prob_wage(workertype,task,task1,phi,salary,V1,V2,V3,V0,V,edu,p,age, year0, year1, year2, year3, year4, year5, &
                      year6, year7, year8, year9, wten)

INTEGER, INTENT(IN):: task, task1, salary,edu,p,age, year0, year1, year2, year3, year4, year5, &
                      year6, year7, year8, year9, workertype 
REAL(dp), INTENT(IN):: phi, V1,V2,V3, V0, V, wten
REAL(dp):: prob_wage
REAL(dp), PARAMETER:: tiny2=1.0d-20
REAL(dp):: arg3(3)
INTEGER:: age1, yearwage
REAL(dp):: constL1, constL2, constL3, slopeL1



!--------------------------------------
!1. Create the variables of interest
!--------------------------------------
age1 = age 

IF(year0 .ne.0) THEN
 yearwage=70

ELSE IF(year1 .ne.0) THEN
 yearwage=71


ELSE IF(year2 .ne.0) THEN
 yearwage=72

ELSE IF(year3 .ne.0) THEN
 yearwage=73

ELSE IF(year4 .ne.0) THEN
 yearwage=74

ELSE IF(year5 .ne.0) THEN
 yearwage=75

ELSE IF(year6 .ne.0) THEN
 yearwage=76

ELSE IF(year7 .ne.0) THEN
 yearwage=77

ELSE IF(year8 .ne.0) THEN
 yearwage=78

ELSE IF(year9 .ne.0) THEN
 yearwage=79

END IF



!=======================
IF(workertype==1) THEN
!=======================
constL1 = s1_5
constL2 = s_002
constL3 = s_003
slopeL1 = d5

!============================
ELSE IF(workertype==2) THEN
!============================
constL1 = s1_52
constL2 = s_0022
constL3 = s_0032
slopeL1 = d5t2

!============================
ELSE IF(workertype==3) THEN
!============================
constL1 = s1_53
constL2 = s_0023
constL3 = s_0033
slopeL1 = d5t3

!============================
ELSE IF(workertype==4) THEN
!============================
constL1 = s1_54
constL2 = s_0024
constL3 = s_0034
slopeL1 = d5t4

!=======
END IF
!=======



IF( wten < 4.0_dp ) THEN


 IF (  (  constL1  + slopeL1*phi +tencoef1*wten + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp &
     +w_74*year4 +w_75*year5 +w_76*year6 +w_77*year7 +w_78*year8  +w_79*year9 ) > tiny2) THEN



 arg3(1) =   constL1 + slopeL1*phi +tencoef1*wten + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp &
     +w_74*year4 +w_75*year5 +w_76*year6 +w_77*year7 +w_78*year8  +w_79*year9 

			  

 ELSE 

 arg3(1) = tiny2


 END IF



IF (  (    constL2 + slopeL1*phi + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp &
      +w2_74*year4 +w2_75*year5 +w2_76*year6 +w2_77*year7 +w2_78*year8  +w2_79*year9 ) > tiny2) THEN


 arg3(2) =  constL2 + slopeL1*phi  + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp  &
      +w2_74*year4 +w2_75*year5 +w2_76*year6  +w2_77*year7 +w2_78*year8  +w2_79*year9  

 
 IF (wten==3.0_dp) THEN

  arg3(2) =  constL2 + slopeL1*phi + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp  &
      +w2_74*year4 +w2_75*year5 +w2_76*year6 +w2_77*year7 +w2_78*year8  +w2_79*year9  

 END IF

 

 ELSE 

 arg3(2) = tiny2


 END IF
 
 
 IF (  (  constL3 + slopeL1*phi + s3_3*REAL(edu,dp)  +s3_1*age1 +s3_2*age1**2.0_dp &
     
      +w3_74*year4 +w3_75*year5 +w3_76*year6 +w3_77*year7 +w3_78*year8  +w3_79*year9  ) > tiny2) THEN


 arg3(3) = constL3 + slopeL1*phi + s3_3*REAL(edu,dp) +s3_1*age1 +s3_2*age1**2.0_dp &
     
      +w3_74*year4 +w3_75*year5 +w3_76*year6 +w3_77*year7 +w3_78*year8  +w3_79*year9  

     	

 ELSE 

 arg3(3) = tiny2


 END IF


ELSE IF( wten >= 4.0_dp ) THEN



 IF (  (  constL1 + slopeL1*phi -tencoef1*wten  + s_3*REAL(edu,dp)+ s_1*age1 + s_2*age1**2.0_dp &
      
	   +  w_74*year4 +w_75*year5 +w_76*year6 +w_77*year7 +w_78*year8  +w_79*year9   ) > tiny2) THEN


 arg3(1) =   constL1 + slopeL1*phi -tencoef1*wten  + s_3*REAL(edu,dp) + s_1*age1 +s_2*age1**2.0_dp &
       
	   +  w_74*year4 +w_75*year5 +w_76*year6 +w_77*year7 +w_78*year8  +w_79*year9 

    
	  

 ELSE 

 arg3(1) = tiny2


 END IF



IF (  (  constL2 + slopeL1*phi + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp  &
       
	   + w2_74*year4 +w2_75*year5 +w2_76*year6 +w2_77*year7 +w2_78*year8  +w2_79*year9 ) > tiny2) THEN


 arg3(2) =  constL2 + slopeL1*phi + s_3*REAL(edu,dp)  +s_1*age1 +s_2*age1**2.0_dp  &
      
	   + w2_74*year4 +w2_75*year5 +w2_76*year6 +w2_77*year7 +w2_78*year8  +w2_79*year9


 
  IF( wten >= 5.0_dp ) THEN	
     
  	arg3(2) =  constL2 + slopeL1*phi + s_3*REAL(edu,dp) +s_1*age1 +s_2*age1**2.0_dp&
	   
	     + w2_74*year4 +w2_75*year5 +w2_76*year6 +w2_77*year7 +w2_78*year8  +w2_79*year9 
  
  	END IF	


 ELSE 

 arg3(2) = tiny2


 END IF
 
 
 
 IF (  ( constL3 + slopeL1*phi + s3_3*REAL(edu,dp) +s3_1*age1 +s3_2*age1**2.0_dp &
      
      +w3_74*year4 +w3_75*year5 +w3_76*year6  +w3_77*year7 +w3_78*year8  +w3_79*year9  ) > tiny2) THEN

 
 
 arg3(3) = constL3 + slopeL1*phi + s3_3*REAL(edu,dp) +s3_1*age1 +s3_2*age1**2.0_dp &
      
      +w3_74*year4 +w3_75*year5 +w3_76*year6 +w3_77*year7 +w3_78*year8  +w3_79*year9  

  IF (wten==6.0_dp) THEN

	  arg3(3) = constL3 + slopeL1*phi + s3_3*REAL(edu,dp) +s3_1*age1 +s3_2*age1**2.0_dp &
      
      +w3_74*year4 +w3_75*year5 +w3_76*year6 +w3_77*year7 +w3_78*year8  +w3_79*year9 
  
  END IF


 ELSE 

 arg3(3) = tiny2


 END IF

END IF


!-------------------------------------------------
!2. Compute the probability of the observed wage
!-------------------------------------------------

!=======================
IF(workertype==1) THEN
!=======================


!+++++++++++++++++++++++++++++++++
IF(task==-2 .or. salary==-2) THEN
!+++++++++++++++++++++++++++++++++
prob_wage = 1.0_dp


!++++++++++++++++++++++
ELSE IF (task==1) THEN 
!++++++++++++++++++++++



prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp *(sigma1_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /( sigma1_u**(2.0_dp)+sigma_eps**(2.0_dp) )  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(1)     )**(2.0_dp))


!++++++++++++++++++++++
ELSE IF (task==2) THEN 
!++++++++++++++++++++++

prob_wage =  (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  (sigma2_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma2_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &

	        * (  DLOG(REAL(salary,dp)) - arg3(2)     )**(2.0_dp))



!+++++++++++++++++++++++
ELSE IF (task>=3) THEN 
!+++++++++++++++++++++++


prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  ( sigma3_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /( sigma3_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(3)     )**(2.0_dp))


ELSE 
PRINT*, "Error in computing the probability of the observed wage for type I"
END IF



!============================
ELSE IF(workertype==2) THEN
!============================


!+++++++++++++++++++++++++++++++++
IF(task==-2 .or. salary==-2) THEN
!+++++++++++++++++++++++++++++++++
prob_wage = 1.0_dp



!++++++++++++++++++++++
ELSE IF (task==1) THEN 
!++++++++++++++++++++++


prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp &
           *  (sigma12_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma12_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(1)    )**(2.0_dp))



!++++++++++++++++++++++
ELSE IF (task==2) THEN 
!++++++++++++++++++++++

prob_wage =  (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  (sigma22_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma22_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &

	        * (  DLOG(REAL(salary,dp)) - arg3(2)   )**(2.0_dp))


!+++++++++++++++++++++++
ELSE IF (task>=3) THEN 
!+++++++++++++++++++++++

prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  ( sigma32_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(  sigma32_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(3) )**(2.0_dp))




ELSE 
PRINT*, "Error in computing the probability of the observed wage type II"
END IF





!===========================
ELSE IF(workertype==3) THEN
!===========================
   
 
!+++++++++++++++++++++++++++++++++
IF(task==-2 .or. salary==-2) THEN
!+++++++++++++++++++++++++++++++++
prob_wage = 1.0_dp



!++++++++++++++++++++++
ELSE IF (task==1) THEN 
!++++++++++++++++++++++


prob_wage = (1.00_dp/( DSQRT((sigma13_u**2.0_dp+sigma_eps**(2.0_dp))*2.0_dp*3.14159265358979323846_dp) ))&
	        *DEXP(   ( -0.5_dp /(sigma13_u**2.0_dp+sigma_eps**(2.0_dp)) ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(1)    )**2.0_dp  )


!++++++++++++++++++++++
ELSE IF (task==2) THEN 
!++++++++++++++++++++++

prob_wage =  (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  (sigma23_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma23_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &

	        * (  DLOG(REAL(salary,dp)) - arg3(2)   )**(2.0_dp))




!+++++++++++++++++++++++
ELSE IF (task>=3) THEN 
!+++++++++++++++++++++++


prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  ( sigma33_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /( sigma33_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(3)  )**(2.0_dp))




ELSE 
PRINT*, "Error in computing the probability of the observed wage for type III"
END IF



!=============================
ELSE IF(workertype==4) THEN
!=============================


!+++++++++++++++++++++++++++++++++
IF(task==-2 .or. salary==-2) THEN
!+++++++++++++++++++++++++++++++++
prob_wage = 1.0_dp



!++++++++++++++++++++++
ELSE IF (task==1) THEN 
!++++++++++++++++++++++


prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp &
            *  (sigma14_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma14_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(1)   )**(2.0_dp)  )


!++++++++++++++++++++++
ELSE IF (task==2) THEN 
!++++++++++++++++++++++

prob_wage =  (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  (sigma24_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /(sigma24_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &

	        * (  DLOG(REAL(salary,dp)) - arg3(2)  )**(2.0_dp))


!+++++++++++++++++++++++
ELSE IF (task>=3) THEN 
!+++++++++++++++++++++++


prob_wage = (1.00_dp/(DSQRT(  2.0_dp*3.14159265358979323846_dp&
            *  ( sigma34_u**(2.0_dp)+sigma_eps**(2.0_dp))  )))&
	        *DEXP(   (-0.5_dp /( sigma34_u**(2.0_dp)+sigma_eps**(2.0_dp))  ) &
	        * (  DLOG(REAL(salary,dp)) - arg3(3)  )**(2.0_dp))




ELSE 
PRINT*, "Error in computing the probability of the observed wage for type IV"
END IF

!======
END IF
!======

END FUNCTION prob_wage

END MODULE mod_wage